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SUMMARY 


A new higher order triangular plate-bending finite element is presented 
in this paper which possesses high accuracy for practical mesh subdivisions 
and which uses only translations and rotations as grid point degrees of free- 
dom. The element has 18 degrees of freedom (d.o. f . ), viz. , the transverse 
displacement and two rotations at the vertices and mid-side grid points a:' the 
triangle. The transverse displacement within the element is approximated by 
a quintic polynomial; the bending strains thus vary cubically within the 
element. Transverse shear flexibility is. taken into account in the stiffness 
formulation. Two examples of static and dynamic analysis are included to 
show the behavior of the element. Excellent accuracy is achieved in all cases. 

This element, designated os TR-18, is demonstrated to be an ideal candi- 
date for generation of a family of plate and shell elements for inclusion into 
NASTRAN. The following elements are specifically mentioned in this context, 
viz., (i) triangular plate element, (ii) quadrilateral plate element, 

(iii) curved triangular shell element, (iv) curved quadrilateral shell element 
and (v) plates with membrane-bending coupling and multilayered plates. The 
present paper describes the detailed theoretical derivations for the afore- 
mentioned elements. In addition, the behavior of the TR-18 element and 
associated quadrilateral plate element is illustrated by two sample problems. 
Comparisons with existing elements in the literature and the present NASTRAN 
quadrilateral elements are shown. 


INTRODUCTION 


NASTRAN presently (Level 15*5) has, in all, a total of nine different 
forms of plate elements in two different shapes (triangular and quadrilateral). 
The present NASTRAN basic bending element, TRBSC, the basic unit from which 
the bending properties of the other plate elements are formed, uses a cubic 
displacement field (with the x^y term omitted). This constrains the normal 
slope (on the exterior edges of the TRPLT bending element) to vary linearly, 
which in turn makes the element overly stiff. A need thus exists for 
a more accurate plate bending element for NASTRAN. 

A brief review of some of the more important plate bending elements is 
now made. Formulations of triangular plate bending finite elements were 
given as long ago as 1966 by Clough and Tocher (ref. l) and by Bazeley et al. 
(ref. 2). The conforming elements presented therein allow only a linear 
variation of slope normal to an edge and have since been found to be overly 
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stiff, whereas the nonconforming element given in ref. 2 uses a cubic poly- 
nomial for transverse displacement and is not of very high accuracy. Improve- 
ments to these elements have been made by using higher degree polynomials for 
transverse displacements; indeed elements of very high accuracy have been 
reported by Argyris (ref. 3), Bell (ref. h ) and Cowper et al. (ref. 5) using 
quintic polynomials for the displacements field. But these elements have 
strains, curvatures and/or higher order derivatives of displacements as grid 
point degrees of freedom (d.o.f.) which lead to an inconsistency when abrupt 
thickness or material property variation occurs. That is to say that the con- 
tinuity of strains and curvatures implied by their use as degrees of freedom 
at grid points is violated wherever concentrated loads, changes in slope, 
changes in thickness, or connections to other structures occur. In short, the 
proper use of elements that assume continuity of strains and curvatures is 
restricted to regions where discontinuities do not occur. Further, the ex- 
istence of higher order derivatives makes it difficult to impose boundary con- 
ditions on these and indeed the simple interpretation of energy derivatives 
as "nodal forces" disappears (ref. 6). Bell has also developed another element 
in ref. k , designated T-15 by him, which has only displacements and rotations 
as degrees of freedom. But it has a major drawback in that not all grid points 
of the element have the same d.o.f.; consequently, it becomes difficult, if 
not impossible, to consider connections of this element with other finite 
elements. Thus the practical use of the T-15 element in general purpose 
programs is severely limited. 

A need still exists to develop a new accurate plate bending finite element 
that has the advantages of the accuracy associated with a high order displace- 
ment polynomial- but does not have the disadvantages discussed above and is 
therefore suitable for inclusion in general purpose computer programs like 
NASTRAN. 

In this paper, a triangular element and an associated quadrilateral 
element are developed that use only displacements and rotations as grid 
point degrees of freedom and use a quintic polynomial for lateral displace- 
ment. The quadrilateral element is formed by four triangular elements. The 
stiffness, consistent mass and load matrices of the separate triangles are 
evaluated and added by the direct stiffness technique to form the respective 
matrices for the quadrilateral. The terms associated with the internal grid 
points are then eliminated by static condensation. None of the elements 
discussed in references 1 to 5 possess the property of transverse shear 
flexibility. This has been taken into account in the present paper by a 
procedure based on that used in NASTRAN (ref. 7).* The components of transvers 
shear strain are quadratic functions of position. Convergence to the limiting 
case of zero transverse shear strain is uniform. 

In addition, three elements, viz., (i) a curved triangular shell element, 
(ii) a curved quadrilateral shell element, and (iii) a multilayered plate 
element can be derived from the TR-18 element. Together with the quadrilateral 
plate element, these elements constitute the TR-18 family of elements. 

•#. - * 1 1 '■ ■ 

A similar procedure for incorporation of transverse shear flexibility 

into a quartic element was communicated to the author by Dr. R. H. MacNeal of 
MacNeal-Schwendler Corporation. 
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LIST OF SYMBOLS 


{a} Column vector of coefficients 

a,b,c Dimensions of triangular element in local co-ordinates (fig. 1) 

a l’ a 2’ a 3 Coefficients of quint ic polynomial 

[B^] , [B^] , [B^] Matrices relating strains and generalized displacements 

[C] Row vector relating transverse displacement to generalized 

displacement 

[D] Matrix relating bending stresses and bending strains 


t 
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w c 
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X,Y,Z 
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\z* ^yz 

b> 
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{<t>} 

v 

P 

X 

[A] 

b) 

U> 


Plate flexural rigidity, Et3/l2(l - v2) 

Elastic modulus 

Matrix relating interior grid point displacement to exterior 
grid point displacements of quad-element 
Matrix relating transverse shear forces and strains 
Stiffness matrix 
Leng-h of side of plate 
Consistent mass matrix 

Vector of bending, a.nd twisting moments per unit length 
Number of elements per side of plate 

Matrix relating grid-point displacement -vector and vector of 
polynomial coefficients 

Augmented matrix of Q and constraint relations 
Matrix relating vector of polynomial coefficients and grid 
point displacement vector 
Kinetic energy 
Transformation matrices 
Thickness of plate 

Matrix of transformation of strain components 
Strain energy 

Vector of transverse shears per unit length 
Lateral displacement 
Central deflection 

Co-ordinate axes in the local system 
Co-ordinate axes in the global system 
Rotation of xz plane at each grid point 
Rotation of yz plane at each grid point 
Transverse shear strains 

Vector of transverse shear strains 

Column vectors of grid point displacement in local or global 
system 

Inclination of material orientation axis to x-axis 
Displacement vector of quadrilateral element 
Poisson's ratio 

Mass density of plate material 

Non-dimensional parameter of eigenvalues, ptu^L^/D 
Direction cosine matrix of quadrilateral median plane 
Circular frequency of plate vibration 
Bending strains 
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Bending strains 


TRIANGULAR PLATE ELEMENT TR-l8 


In this section of the paper, the derivation of the stiffness matrix, 
consistent load vector and consistent mass matrix of the triangular plate 
element is given. The procedure for the derivation is described in detail in 
reference 3, and hence only essential details arc presented here. 

The element has 18 d.o.f., the transverse displacement and 2 rotations 
at each vertex and at the mid-point of each side. Three additional conditions 
are introduced, viz., the slope normal to each edge (hereinafter called normal 
slope) varies cubically along each edge. This establishes 3 constraint equation 
between the coefficients of the polynomial for displacements, which, together wi 
the 18 d.o.f., uniquely determine the 21 coefficients in the quintic polynomial. 
The variation of deflection along any edge is a quintic polynomial in the 
edgewise co-ordinate; the six coefficients of this polynomial are uniquely 
determined by deflection and edgewise slope at the 3 grid points of the edge. 
Displacements are thus continuous between two elements that have a common edge. 
The normal slope along each edge is constrained to vary cubically; however, 
since the norma.*, slopes are defined only at 3 points along an edge, there is 
no normal slope continuity between 2 elements that have a common edge. The 
element thus belongs to the class of non-conforming elements. The development 
of this element follows closely that of Cowper et al. (ref. 5). 


Element Geometry 

Rectangular cartesian co-ordinates are used in the formulation. An 
arbitrary triangular element Is shown in figure 1, where X, Y, and Z are a 
system of global co-ordinates and x, y, z are the system of local co- 
ordinates for the triangular element. The grid points of the element are 
numbered in counterclockwise direction as shown. The following relationship 
between the dimensions of the triangular element a, b, c, the inclination 
0 between the X and x axes and the co-ordinates of the vertices of the 
element can be easily derived (see fig. 1) : 


where 


cos 0 * (X^ - X^)/r si* 1 6 * (¥3 " Y^)/r 

2 2 1/2 
r = l(X 3 - X x r ♦ (Y 3 - Y^ ) j 


( 1 ) 

( 2 ) 


a * (X 3 - X^) cos 0 - (Y,. - Y 3 ) sin 0 


{ (X 3 - X 5 )(X 3 - X^ + (Y 3 - * 5 )(Y 3 - Y 1 )}/r 


(3) 
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Displacement Function 

The deflection w(x, ” ) within the triangular element is assumed to vary 
as a quintic polynomial in the local co-ordinates, i.e.. 
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There are 21 constants, a^ to a^^. These are evaluated as follows: 

The e.lement has 18 d.o.f. At each grid point there are 3 displacement 
components as d.o.f., viz., w, displacement in z-direction, o, rotation 
about the x-axis and 8, rotation about y-axis. The rotations a and 8 
are obtained from the definitions of transverse shear strains y and 


Y , i.e. , 
yz * 


xz 


Y ■ ip ♦ 8 
xz dx 

dv 

Y * t— - a 

yz 3y 


(7) 


It can be shown (ref. 8) that y and y . , and hence a and 8, at any 

xz yz 

grid point can be expresred in terms of the constants a^ to a^. Thus 18 

relations between grid point displacement values and the constants are obtained. 
Three constraints among the coefficients in the above polynomial (eq. (6)) are 
now introduced so that the normal slope varies cubically along each edge. It 
is clear that the thr'e constraint equations will involve only the coefficients 
of the fifth degree terms in equa.xon (6), since the lower degree terras 
satisfy the condition of cubic normal slope automatically. Moreover the con- 
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ditions depend only on the direction of nn Mge and not on its position, Along 
the edge defined by grid points 1 and 3, where y = 0, the condition of cubic 
normal slope requires that 


a 17 = 0 (8) 

It can be shown (ref. 8) that the condition for cubic variation of normal slope 
alone edge 1-5 is 

5b^c + (Ub 3 c 2 - b'ta^ + (3b 2 c 3 - 2b^c)a 1 g + (2bc^ - 3b^c 2 )a 1 ^ + 

(c 5 - l*b 2 c 3 )a 20 - 5bc a 21 = 0 (9) 


and the condition for cubic variation of the normal slope along the edg° 3-5 
(see fig. 1) is 

5a^c a^g + (-4a 3 c 2 * a^a^ + (3a^c 3 - 2a**c)a 1 g + (-2 ac** + 3& 3 c‘ L )a 1 g + 

(c 5 - 4a 2 c 3 ) a 2Q + 5ac^a 21 = 0 (10) 


The 18 relations between grid point displacements (w, a and 8 at each 
of the six grid points) and the coefficients of the polynomial, together with 
the three constraint equations (8), (9), and (10 >, uniquely determine the 
coefficients a^ to a^. The following equations can therefore be written: 


{6} - IQ] {a} 


( 11 ) 


and 


{a} = [SJ {5} 


( 12 ) 


where IQ] is the 18 x 21 matrix involving the co-ordinates of grid points 
substituted into the function w (eq. (6)) and the appropriate expressions 
of a and 8; {a} is the column vector of coefficients a^ to a 21 , and [S] 

is a 21 x 18 matrix and consists of the first 18 columns of the inverse of an 
augmented matrix o'* [Q] and the three constraint equations (o), (9)» and (10). 


Stiffness Matrix 

The following relationships are obtained from the theory of deformation 
for plates (ref. 9) . In the present notation, the curvatures are defined by 
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lending and twisting moments are related to curvatures by 


( 13 ) 
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(14) 


where [D] is, in general, a full symmetric matrix of elastic coefficients. 

Shear forces (and hence shear strains) are proportional to the third 
erivatives of the displacements. Since the displacement within the element 
's assumed to vary as a quintic polynomial, shear strains are expressed by a 
jjadratic polynomial as follows: 

Y x * b i + V + b^y + b^x* + b^xy + b g y 2 ^ 


2 2 

Y y = C 1 + C 2 X + c 3 y + C 4 X + c 5 xy + c 6 y 




(15) 


The shear forces 1 , V are related to y , v by 

x y 'x* 'y J 


= t* [G]< 




( 16 ) 


'here G is in general a full 2x2 symmetric matrix and t* is an effective 
thickness of the element. 


It can be shown that b, to b^ and c, to c, can be expressed in 

1 6 i 6 r * 

J xz lean be 

IVJ 


terms of the coefficients a^ to (ref. 8) and hence 
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expressed as 


{y} = [B-jJ {a} 


(IT) 


where [B^] is as given in reference 8. The curvature (x) is now split into 
2 parts, i.e. , 


where 
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(19) 


It follows that (x^l is the vector of curvatures in the absence of trans- 
verse shear and (x^) is the contribution of transverse shear to the vector of 
curvatures. Now ,X-^) and {x 2 } are expressed in terms of generalized co-ordinate 
{a} as 


* X 1* = 


( 20 ) 


and 


(x 2 } = [B 3 ] (a) 


where [B^] and [B^] are given in reference 8. Thus, 

M = {X x } + (X 2 ) = ([B 2 ] + [B 3 1) {a} 

The generalized stiffr.ess matrix can be obtained as (ref. 8) : 
lK, s*n ■ II {( 1E 9 4 ‘Vf'^nv * '“s’) 4 IB / dx dy 


( 21 ) 


( 22 ) 


(23) 


462 


. • : «.» • yV f, 


t- 

't 


The 


element stiffness matrix in the local co-ordinate system, 


lyirtue of equation (12), 
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IKK = IS) 1 [K) sen IS] 
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f The element stiffness matrix in the global co-ordinate system, [K] , is 


[K] = (T 2 r [K] e [T 2 1 


(25) 


■tfr 


i 


where [T ] is the transformation matrix of displacement vectors from global 
to local co-ordinates of element. 

The evaluation of the elements of the generalized stiffness matrix, [K] gen 

of equation (23), in closed form is, though straightforward, very tedious. 

This is due to the lengthy expressions involved in the triple matrix products. 
The integration involved in equation (23) is now split into 5 integrals as 
follows : 

[K] _ = // IBJ T [D] [B 1 dx dy 


■gen 


+ // [B ] T [D] [B ] dx dy + ff [Bg] 1 [D] [Bg] dx dy 


+ ff (B ] T [D] [B ] dx dy + // [B ± ] l G l 1^1 dx ^ 


(26) 


The first term ff [B 2 ] T [D] [Bg] dx dy is evaluated in closed form; the 

other four terms are evaluated using numerical integration. The numerical in- 
tegration formulas used are listed in ref. 8. If the plate is assumed o be 
rigid in transverse shear, the matrices [B^] and [B^] are null and the las 
four terms of equation ( 26 ) vanish. 

Consistent Mass Matrix 

It can be shown that the generalized consistent mass matrix is (ref. 8) 


[M] = pt ff [C] [C] dx dy 

gen 


(27) 
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2 p 5 

where [C] = (1 x y x xy y y ]. 

The mass matrix can be transformed to element co-ordinates and global 
co-ordinates by the same transformations as those used for stiffness matrix. 
Thus , 


[M] 0 = [S P [M] gen [S] (28) 

and 

[Ml = [T ] T [Ml [T ] (29) 

g 2 e 2 

where the subscripts e and g on [M] stand for element and global system, 
respectively. 


Consistent Load Vector 

It can be shown that the generalized consistent load vector is, (ref. 8) 

[p] „ = ff t c ^ T dx d y (30) 

gen 

where q is the distributed loading. 

The consistent load vect^ :.n now be transformed to element and global 
co-ordinates by 


[P] _ = [S1 T [PI 

(31) 

e 

gen 


w, - 1 

[t/ [Pl e 

(32) 


THE QUADRILATERAL PLATE ELEMENT 


The quadrilateral element is formed from four of the triangular elements 
just described. Two arrangements of the quadrilateral element are shown in 
Figures 3(a) and 3(b). 

The quadrilateral element has eight grid points on its edges. In the 
arrangement of the quadrilateral element shown in Figure 3(a), which will be 
designated as QUAD1, the quadrilateral is divided, first into 2 triangles by 
one diagonal and then again into 2 more triangles by the other diagonal. In 


U6'i 


each case one additional grid point, at the mid-point of the diagonal, is intro- 
duced j the stiffness, mass and load matrices of the triangles are evaluated and 
added and the terms associated with the internal grid point are eliminated by- 
static condensation. The stiffness, mass and load matrices of the quadrilateral 
element are obtained by adding one-half the contribution of each case. In the 
arrangement of the quadrilateral element shown in Fig. 3(b), designated as 
QUAD5, five additional grid points are introduced internally so that the quad- 
rilateral is divided into four triangular elements. The eight grid noints on 
the edges are numbered 1 to 8. Grid point 9 is located at the intersection 
of lines joining mid-points of opposite edges. Grid points 10 to 13 are 
located at the middle of the lines joining grid point 9 to each of the corners 
of the quadrilateral. The stiffness, mass and load matrices of the triangular 
elements are evaluated, as described previously, and added by the direct stiff- 
ness technique to form the respective matrices for the quadrilateral. The 
internal grid points are then eliminated by static condensation. 

In a preliminary operation, the grid points of the quadrilateral are ad- 
justed to lie in a median plane. The median plane is selected to be parallel 
to, and midway between, the diagonals of the quadrilateral. The adjusted quad- 
rilateral is the normal projection of the given quadrilateral on the median 
plane. The short line segments joining the corners of the original and pro- 
jected quadrilateral elements are assumed to be rigid in bending and extension. 
The quadrilateral element and its projection onto the median plane is shown in 
Fig. 3(c). 


FORMULATION AND SOLUTION OF EQUATIONS 


The global stiffness matrices, load vectors, and mass matrices for the 
complete structure modeled by these elements are assembled from the correspond- 
ing matrices of the individual elements by standard methods (ref. 6) to 
form the matrix equation 


[K] {U} = (P) 


(33) 


Because the d.o.f. at grid points consist of displacements and rotations, 
it presents no difficulty to specify the appropriate geometric boundary 
conditions at any irregular and/or complex boundary. After the boundary 
conditions are applied, the matrix equation (33) is solved by Gaussian 
elimination to obtain the global displacement vector {Uh 


DISCUSSION OF RESULTS 


The triangular and quadrilateral elements are used to solve two problems 
in statics and dynamics of thin isotropic plates. Only the results for the 
simply supported plate are presented here; the interested reader may consult 
ref. 8 for details. The problem analyzed is that of the statics and dynamics 
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of e square plate with edges simply supported. All calculations were carried 
out on the CPC 6400/6600 series of computers with SCOPE operating system of 
the Langley Research Center. Single precision arithmetic was used xhroughout. 

A value of Poisson’s ratio of 0.0 is used in all problems. It is mentioned 

in this context that other finite-element analyses in the literature use 0.3 ; 

as the value of Poisson’s ratio. i 
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Static Analysis of a Square Plate 

The arrangement of the finite elements in a quarter of the square plate 
is shown in Fig. 4. The number of subdivisions of the edge of the square is 
denoted by N. Due to symmetry, only one-quarter of the plate is analyzed. 

The calculated values of the deflection at the center of the simply supported 
plate are given in Table 1 and compared with the exact solution given by 
Timoshenko (ref. 9)* These values together with other known finite element 
analyses available in the literature (refs. 3> 4, 5 and 10) are also compared 
in Figures 5 and 6 in plots of deflection versus mesh size using a linear scale 
for N~^. 

As seen from table 1, the "Q" arrangement is found to give better results 
than the ”P” arrangement for the uniformly distributed loading; however, the 
"P ff arrangement is found to be better, in general, for concentrated loads. For 
the clamped plate, the ff P" arrangements are found to be slightly better than 
the "Q” arrangements, as noted from ref. 8. For the quadrilateral element, 
QUAD1 is found to be superior bo QUAD 5 . 

The high accuracy achieved with the present elements (triangular and 
quadrilateral), even for the coarsest mesh, is evident from Table 1 and 
Figures 5 and 6 for the simply supported plate. In the case of the clamped 
plate, the results for the coarsest grid are not as accurate as in the case of 
the simply supported plate (ref. 8); however, as the element size is decreased 
the values of deflection obtained with the present elements approach very 
rapidly the exact results. 


Free Vibration of a Square Plate 

The natural frequencies of a simply supported square plate wei\ letermined 
using the triangular and quadrilateral elements. The non-dimensional eigen- 
values are 

A = ptccPlA/D (3*0 


p = mass density 
t = thickness of plate 
(*) = circular frequency 
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L = length of side of square plate 

D = Et3/l2(l - v ^) 9 the flexural rigidity of the plate. 

The exact eigenvalues for the simply supported plate are given by 

X = {r 2 + s 2 ) 2 it 1 * (35) 


where r and s refer to the number of half-waves parallel to the edge directions. 

The lowest 6 values obtained using the present elements arid the exact 
results are shown in Table 2. The eigenvalue problems were solved using a 
Jacobi routine that produced the complete set of eigenvalues and eigenvectors. 
Consistent mass matrix was used for treatment of inertia. It is seen that the 
lowest eigenvalue is calculated to within 1 % of exact result. Good agreement 
is noticed for higher eigenvalues as well. 


THE TR-18 FAMILY OF ELEMENTS 


A number of finite element formulations for doubly curved shells are 
presently available, the notable among them being the works of Ahmad, Irons 
and Zienkiewicz (Ref. 11), bonnes, Dhatt, Giroux, and Robichand (Ref. 12), 
Strickland and Loden (Ref. 13), Key and Beisinger (Ref. 1*0, Dhatt (Ref. 15 )» 
and Olson and Lindberg (Ref. l6). Some of these have neglected transverse 
shear deformations whereas some others use sub-triangles and/or second and 
higher order derivatives of the displacements of the element as degrees of 
freedom, thus complicating the formulation. A need still exists for an 
accurate shell element that has only translations and rotations as d.o.f. 

Such shell elements can be derived using the TR-18 plate element; the 
formulation presented here is simple and includes transverse shear deforma- 
tions; it is based on the linear shear deformation theory of thin shells as 
given by Washizu (Ref. 17). 

Using shallow shell theory, flat plate elements can be easily converted 
to curved shell elements. The linear strain triangular membrane element, 
known as TRIM6 in the literature, can be combined with the TR-l8 plate element 
to develop a doubly curved shallow shell triangular element. The surface of 
the shell will be approximated by a quadratic polynomial of the position 
coordinates of the base triangle. By a procedure analogous to that discussed 
for the quadrilateral plate element, a quadrilateral shallow shell element 
can be developed. Multilayered plates, and plates with coupled membrane and 
bending deformations, can be designed using TR-18 plate elements. 
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Curved Triangular Shell Element 


Fig. 7 shows a differential element dA on the middle surface of the 
doubly curved shell with orthogonal curvilinear surface co-ordinates E^, 

A right handed cartesian co-or be system X, Y, Z is also shown. 

In Fig. 8 and Fig. 9 the curved triangular shell element is shown in basic and 
local coordinate systems. The differential surface element is expressed as 


dA = a l a 2 d ^l d ^2 


(36) 


where and are the Lame parameters. 


If the surface z(x,y) of an element is shallow, the following relations 
are valid 
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The set of orthogonal curvilinear co-ordinates (E^, E^» o ver the surface 

of the shallow element dA can be replaced by a set of shallow cartesian 
co-ordinates (x,y,z) where 




= x C 2 = y 

(38) 

and Lam4 

parameters 

a i = a 2 ^ 1 

(39) 

From eq. 

(36), (38) 

and (39). dA = dx dy 

(40) 

The 

curvatures 

of the shallow element can then be approximated by 
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is the vector from the origin of thi basic co-ordinate system to 


the origin of the local co-ordinate system. The distances X , Y , Z are 

'JOG 

not involved in the calculation of the stiffness matrix of the element since 
only the differences of co-ordinates are used; hence they are discarded. The 
inversion of equation (44) yields 
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A _ 

11 
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It may be seen that A.^, A,^ 


and A , 


are components in X, Y and Z 


directions of a unit vector along the x direction; and so on for A^, A , e ^ Ci 

An analytical description of the surface of the element suitable for 
application of shallow shell theory is obtained by assuming that the elevation 
of the shell middle surface may be expressed as a quadratic polynomial in the 
local co-ordinates of the element, i.e., 

z(x,y) = f 1 + f 2 x + f^y + f^x 2 + f ? xy + f^y 2 (46) 
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This implies that the. shell element has constant curvatures and is consistent 
with the approximations of shallow shell theory. Knowing the co-ordinates 
x, y, z of the six points of the triangular 'lement, the constants f to 

fg can he evaluated. 

Symbolically 


{ z } = [1^] {f} 

(47) 

U‘} = [Q 1 ]” 1 (z) 

(48) 


where [Q^] is a 6 x 6 matrix of the co-ordinates of the six points of the 
element substituted into equation (46). 


D egrees of freedom and assumed displacement function .- The element has 30 
degrees of freedom (d.o.f.), with 5 d.o.f. per grid point. These are the 
three translations u. v, w in the x, y, and z directions and the rota- 
tions of the xz and yz planes, a and 3. The displacements u, v, w 
are positive in the positive co-ordinate directions j the slopes are positive 
when they cause compression at the top of the surface. The u and v d.o.f. 
are assumed to vary ever the element by a full quadratic polynomial of local 
co-ordinate", as follows: 


2 2 

u = a^ + a ? x + a^y + a^x + a^xy + a^y 


v = a ? + a 8 x + a ? y + a 1Q x + a„xy + a 10 y 


11 12 


(49) 

( 50 ) 


The deflection w will be defined by a quintic polynomial as in equation (6). 
The coefficients a! to of equation (6) will be renumbered a-^ to a^ 

respectively. The 33 coefficients a^ to a^ can be uniquely determined 
from the 30 d.o.f. of the shell element (5 d.o.f. each at six grid points) 
together with the 3 constraint equaticus (8) , (9) , and (10) . 


Strain-displacement relations .- The expressions for transverse shear 
strains and bending strains for the curved shell element are the same as those 
for the TR-18 element (eqs. (7) and (13)). The membrane strains are 
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>r^; 


e 

x 


e 

y 


_ 3u 3 ? z 

' ax “ W a 2 
3x 


av a 2 z 





_ 3u 3v ? 3 2 z 

e xy 3y 3x ” w 3x3y 


( 51 ) 


Stiffness matrix .- The stiffness matrix can be evaluated by the standard 
procedures (ref. 6). The element can then be tested against other elements 
(refs. 11 to 16) for suitability as well as accuracy. At the time of writing 
of this paper, the calculations for the element have not been completed. 


Curved Quadrilateral Shell Element 


A curved quadrilateral shell element can be constructed from the curved 
triangular shell elements by a procedure analogous to that of the construction 
of the quadrilateral plate element from the TR-18 element. 


Plates With Membrane-Bending Coupling 

Plates with coupled membrane and bending deformations and multilayered 
plates (fig. 10) can be analyzed by means of the elements presented earlier 
herein. Multilayered plates will produce coupling between membrane and 
bending deformations when the plate is not symmetrical with respect to its 
middle surface. A general form of the coupled stress-strain relationship 
can be expressed as 

r 


(F) 


[A] 

13] 

0 


U > 

in 

UM} 


IB) T 

[D] 

0 

< 

<X> 

(V) 

J 


0 

— 

0 

t*[0] 


(y) 

L J 


(5&) 


where 


{F} is a vector of membrane force components 


F , F 
y xy 

(M) is a vector of bending and twisting moments M , M , M 

X Jf 4 VjT 

{V} is a vector of transverse shear components V , V 

* V 
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(e } is a vector of membrane strain components e , e , e 
m x y xy 

{\} is a vector of curvatures x » X * X 

x y xy 

{ Y } is a vector of average transverse shear strain y , y 

x y 

N 

[A] is a 3 x 5 matrix, ^ [0 ] 

k=l 
N 

[21 is a 3 x 3 matrix, £ [G ] 

k=l 
N 

[D] is a 3 x 3 matrix, } [G 1 

L e 
k=l 

[G] is a 2 x 2 transverse shear matrix 

[G 1 is a 3 x 3 matrix of elastic coefficients 
e 

t is the distance to the outer edge of plate (or layer in a multilayer 
plate) from reference surface 

^ the distance to the inner edge of plate (or layer in a multilayered 
plate) from reference surface 

t* is an effective thickness for the element 
The inplane strain vector at any point is 

(e) = {e m ) - z { X } (53) 

where z is the distance from the reference surface. The strain energy of the 
plate element is 

U = | / [{F} T {e m > + {M} T { X } + {V} T {y}]dA (54) 

where the integration is carried out over the surface of the element. The 
stiffness matrix for the triangular and quadrilateral elements can be 
evaluated by the usual procedures (refs. 6 and 8). 


<V - Vi> 


2 2 
t - t 
k k-1 


t k ■ i 
k k-1 


CONCLUDING REMARKS 


New triangular elements and associated quadrilateral elements for plate 
and shell analysis having only displacement and rotations as grid point 
degrees of freedom are described in this paper. The examples presented for 
plate elements demonstrate that high accuracy is achievable using these 
elements for practical subdivisions. 
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f 





The effect of transverse shear deformations is included in the element 
formulation. Transverse shear strains vary quadrat ically within the element; 
convergence to the limiting case of zero transverse shear strain is uniform. 
The present elements are expected to give better approximations than most 
displacement model plate bending elements for solving problems where trans- 
verse sh^ar effects are significant. 

Finally, it is remarked that these elements are ideally suited for 
inclusion into general purpose computer programs due to (i) simplicity of 
formulation, (ii) use of only displacements and rotations as grid point 
degrees of freedom, (iii) high accuracy for practical mesh subdivisions and 
(iv) inclusion of transverse shear flexibility in the element properties. 
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Notation 
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T (T-X8) 
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T (T-21) 
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C-N 

T ( Q-mesh) 

5 

C-P 

T (P-mesh) 

5 

NO 

Q 

7 

CFQ 

Q 

9 

T - Triangular; 0 - Quadrilateral; 
R - Rectangular 
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Figure 5.- Simply supported square plate: central deflection w c 

under central point load P. 
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Figure 6.- Simply supported square plate: central deflection w„ 

under uniformly distributed load q Q , 
















